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Abstract 



(-*. We present an innovative numerical discretization of the equations of inviscid 

,_^ potential flow for the simulation of three-dimensional, unsteady, and nonlinear 

water waves generated by a ship hull advancing in water. 

1 j ' The equations of motion are written in a semi-Lagrangian framework, and 

> the resulting integro-differential equations are discrctizcd in space via an adap- 

^_< tive iso-parametric collocation Boundary Element Method, and in time via im- 

(-H plicit Backward Differentiation Formulas (BDF) with adaptive step size and 

variable order. 

When the velocity of the advancing ship hull is non-negligible, the semi- 
i C i Lagrangian formulation (also known as Arbitrary Lagrangian Eulerian formula- 

tion, or ALE) of the free-surface equations contains dominant transport terms 
CN| which are stabilized with a Streamwise Upwind Pctrov-Galerkin (SUPG) method. 

The SUPG stabilization allows automatic and robust adaptation of the spa- 
tial discretization with unstructured quadrilateral grids. Preliminary results are 
_~. presented where we compare our numerical model with experimental results on 

(f) a Wigley hull advancing in calm water with fixed sink and trim. 

C*~) Keywords: unsteady ship-wave interaction; nonlinear free-surface problems; 



semi-Lagrangian formulation; arbitrary Lagrangian Eulerian formulation; 
boundary clement method 



1. Introduction 



5— i 

C^ The use of computational tools to predict hydrodynamic performances of 

ships has gained a lot of popularity in recent years. Models based on potential 

flow theory have historically been among the most successful to assess the wave 

pattern around a ship hull in presence of a forward ship motion. 

In this framework, the assumptions of irrotational flow and inviscid fluid 

reduce the Navier Stokes incompressibility constraint and momentum balance 
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equations to the Laplace's and Bernoulli's equations, denned on a moving and 
a-priori unknown domain. 

This boundary value problem is tackled with a Mixed Eulerian-Lagrangian 
approach, in which the Eulcrian field equations are solved to obtain the fluid 
velocities, which are then used to displace in a Lagrangian way the free-surface, 
and update the corresponding potential field values [24]. In this framework, 
the Eulerian problem is expressed in boundary integral form, and it is typically 
discretized using the Boundary Element Method (BEM) . The velocity field and 
Bernoulli's equation provide a kinematic boundary condition for the Lagrangian 
evolution of the free-surface, and a dynamic boundary condition for the evolu- 
tion of the potential field. 

Numerical treatment of the Lagrangian step usually relies on accurate recon- 
structions of the position vector and of potential field gradients. In presence of a 
forward motion, these reconstructions may lead to instability in the time advanc- 
ing scheme for the free-surface discretization, as well as for the corresponding 
potential field values. A smoothing technique is often adopted to reduce the 
sensitivity of the discretization on the reconstruction of the full velocity field, 
at the cost of introducing an artificial viscosity in the system. 

An alternative cure was presented by Grilli et al. [13] who developed a high 
order iso-parametric BEM discretization of a Numerical Wave Tank to simulate 
overturning waves up to the breaking point on arbitrarily shaped bottoms. The 
use of a high order discretization bypasses the problem of reconstructing the 
gradients, and is very reliable when the numerical evolution of the free-surface 
is done in a purely Lagrangian way. 

Ship hydrodynamics simulations, however, are typically carried out in a 
frame of reference attached to the boat, requiring the presence of a water current 
in the simulations. In this case a fully Lagrangian approach leads to downstream 
transportation of the free-surface nodes, or to their clustering around stagnation 
points, ultimately resulting in blowup of the simulations (see, for example, [22]). 

Beck [4] proposed alternative semi- Lagrangian free surface boundary condi- 
tions, under the assumption that the surface elevation function is single- valued. 
Employing such conditions, it is possible to prescribe arbitrarily the horizontal 
velocity of the free surface nodes, while preserving the physical shape of the 
free-surface itself. 

However, if there is a significant difference between the water current speed 
and the horizontal nodes speed, stability issues may arise [22], [21]. For this 
reason, the semi-Lagrangian free-surface boundary conditions proposed by Beck 
have been in most cases employed imposing nodes longitudinal speeds equal to 
the water current ones [22], thus reducing the instability at the cost of more 
complex algorithms for mesh management. 

More recently, Sung and Grilli [23] applied an alternative method, combin- 
ing semi-Lagrangian and Lagrangian free surface boundary conditions to the 
problem of a pressure perturbation moving on the water surface. The semi- 
Lagrangian free-surface boundary conditions were used also in the work of Kjcll- 
berg, Contento and Jansson [18], where free-surface instabilities are avoided by 
choosing an earth-fixed reference frame, in which no current speed is needed. 



The drawback of this choice is that in such frame the ship moves with a specified 
horizontal speed, and the computational grid needs to be constantly regenerated 
to cover the region surrounding the hull with an adequate number of cells. 

The purpose of this work is to present a new stabilized semi-Lagrangian 
potential model for the simulation of three dimensional unsteady nonlinear wa- 
ter waves generated by a ship hull advancing in calm water. The resulting 
integro-differential boundary value problem is discretized to a system of nonlin- 
ear differential-algebraic equations, in which the unknowns are the positions of 
the nodes of the computational grid, along with the corresponding values of the 
potential and of its normal derivative. 

Time advancing of the nonlinear differential-algebraic system is performed 
using implicit Backward Differentiation Formulas (BDF) with variable step size 
and variable order, implemented in the framework of the open source library 
SUNDIALS [15]. The collocated and iso-parametric BEM discretization of the 
Laplace's equation has been implemented using the open source C-\ — \- library 
deal.II [3]. 

The computational grid, composed by quadrilateral cells of arbitrary order, 
is adapted in a geometrically consistent way (see [6]) via an a-posteriori error 
estimator based on the jump of the solution gradient along the cell internal 
boundaries. Even when low order boundary elements are used, accurate esti- 
mations of the position vector and potential gradients on the free-surface are 
recovered by a Streamwise Upwind Petrov-Galerkin (SUPG) projection, which 
is used to stabilize the transport dominated terms. The SUPG projection is 
strongly consistent, and does not introduce numerical dissipation in the equa- 
tions. This allows the use of robust unstructured grids, which can be generated 
and managed on arbitrary hull geometries in a relatively simple way. 

The test case considered in this paper is that of a Wigley hull advancing at 
constant speed in calm water. The simulations have been performed using bilin- 
ear elements. The arbitrary horizontal velocity specified in the semi-Lagrangian 
free-surface boundary condition is chosen so the nodes have null longitudinal 
velocity with respect to the hull. The numerical results obtained imposing six 
different Froude numbers are finally compared with experimental results re- 
ported in [20], to assess the accuracy of the model. 

2. Three-dimensional potential model 

The equations of motion that describe the velocity and pressure fields v 
and p of a fluid region around a ship hull arc the incompressible Navier-Stokcs 
equations. Such equations, written in the moving domain fl(t) C K 3 (a simply 
connected region of water surrounding the ship hull) read: 
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where p is the (constant) fluid density, b are external body forces (typically 
gravity and possibly other inertial forces due to a non inertial movement of 
the reference frame), <x = — pi + /i(Vt) + Vt> T ) is the stress tensor for an 
incompressible Newtonian fluid, T(t) := dCl(t) is the boundary of the region of 
interest, and n is the outer normal to the boundary T(t). On the free-surface 
T w (t), the air is assumed to exert a constant atmospheric pressure p a on the 
underlying water, and we neglect shear stresses due to the wind. On the other 
boundaries of the domain, the prescribed velocity v g is either equal to the 
ship hull velocity, or to a given velocity field (for inflow and outflow boundary 
conditions far away from the ship hull itself). 

Equation (la) is usually referred to as the momentum balance equation, 
while (lb) is referred to as the incompressibility constraint, or continuity equa- 
tion. 

In the flow field past a slender ship hull, vorticity is confined to the boundary 
layer region and to a thin wake following the boat: in these conditions, the 
assumption of irrotational and non viscous flow is fairly accurate. The neglected 
viscous effects can be recovered, if necessary, by other means such as empirical 
algebraic formulas, or — better — by the interface with a suitable boundary layer 
model. 

For an irrotational and invishid flow, the velocity field v admits a represen- 
tation through a scalar potential function $(x,i), namely 

v = V$ in n(t). (2) 

In this case, the equations of motion simplify to the unsteady Bernoulli equation 
and to the Laplace equation for the flow potential: 

d ^ + \\V<$>\ 2 + V -^^+(3-x = C{t) infi(t) (3a) 

at 2 p 

A$ = in Q(t) (3b) 

where C(t) is an arbitrary function of time, and we have assumed that all body 
forces could be expressed as b = — V(/3 • aj), i.e., they are all of potential type. 
This is true for gravitational body forces and for inertial body forces due to 
uniform accelerations along fixed directions of the frame of reference. 

In this framework, the unknowns of the problem $ and p are uncoupled, 
and it is possible to recover the pressure by postprocessing the solution of the 
Poisson problem (3b) via Bernoulli's Equation (3a). 



Since arbitrary uniform variations of the potential field produce the same 
velocity field (i.e., V($(aj,t) + C(t)) = V$(x,t)), we can set C(t) = 0, and 
we can combine Bernoulli equation (3a) and the dynamic boundary condition 
on the free-surface (lc) to obtain an evolution equation for the potential field 
$(x,t) on the free-surface T w (t). 

The shape of the water domain Ct(t) is time-dependent and it is part of the 
unknowns of the problem. The free-surface T w (t) should move following the 
fluid velocity v. In ship hydrodynamics, however, it is desirable to maintain 
the frame of reference attached to the boat, and study the problem in a domain 
which is neither fixed nor a material subdomain. Indeed, its evolution is not 
governed solely by the fluid motion, but also by the motion of the reference 
frame and by the motion of the stream of water. In addition, it has to comply 
to the free-surface boundary T w (t) which is the result of the dynamic boundary 
condition (lc). 

A possible solution is to introduce an intermediate frame of reference, called 
Arbitrary Lagrangian Eulcrian (ALE) (see, for example, [10] and Appendix A). 
In the context of potential free-surface flows, this approach is also known as the 
semi- Lagrangian approach [5, 22]. 

2.1. Perturbation potential and boundary conditions 




Figure 1: The domain f!(i) and the different regions in which its boundary dQ(t) is split. 
The grey area behind the hull T h (t), is the absorbing beach portion of the free-surface region 

r w (t). 



We solve Problem (3) on the region Q(t) around the boat, and we move the 
local frame of reference with horizontal velocity Vf(t) which coincides with the 



horizontal boat velocity (which we assume as known) . An additional horizontal 
water stream velocity V s (t), expressed in a global (earth- fixed) reference frame 
is added to the problem. 

Uniform accelerations of the reference frame a,{ can be taken into account 
by incorporating them in the body force term 

j3 ■ x = at ■ x + gz, (4) 

and it is convenient to split the potential $ into the sum between a mean flow 
potential (due to the stream velocity and to the frame of reference velocity) 
and the so called perturbation potential <f> due to the presence of the ship hull, 
namely 



*(x, t) = (V s (t) - V t (t)) ■ x + 4>(x, t) (5a) 

v(x, t) = V$(x, t) = V s {t) - V f (t) + V<j>{x, t) (5b) 

^{x,t) = (a s (t) - af(t)) • x + ^(x,t). (5c) 

The perturbation potential still satisfies Poisson equation 

A0 = O in fl(t), (6) 

and in practice it is convenient to solve for (f>, and obtain the total potential $ 
from equation (5). 

In Figure 1 we present a sketch of the domain fi(£), with the explicit splitting 
of the various parts of the boundary T(t). On the boat hull surface, the non- 
penetration condition takes the form 

^:=V^-n = n-(V b -y s ) on r l (i), (7) 

when expressed in terms of the perturbation potential and Vb is the (given) 
boat velocity. 

On the — horizontal — bottom of the water basin T b , the non penetration 
condition is also applied, namely 

<f>n = on T b (t). (8) 

A possible condition for the — vertical — far field boundary is the homogeneous 
Neumann condition 

0„ = O on T"{t). (9) 

The most remarkable limit of such condition is the fact that it reflects wave 
energy back in the computational domain, thus limiting the simulation time. 
Different boundary conditions can be used to let the wave energy flow outside 
the domain. In Appendix B, we explain in detail the procedure used in our 
computations. 



Applying the potential splitting (5) to Problem (3), and summarizing all 
boundary conditions, we obtain the perturbation potential formulation 

A(f> = in fl(t) (10a) 

d ±+vV<P = -gz + a s -x+ 1 -\Vcf ) \ 2 -^ n ouT w {t) (10b) 

4>(x,0) = M X ) onr(O) (10c) 

0n = n ■ (V b - V„) on T h (t) (lOd) 

n = o onr 6 (t)ur^(t) (10c) 

— Qf 1 = v onr-(t), (lOf) 

where t> = V s — V{ + V</>. Equation (lOf) is a kinematic boundary condi- 
tion and indicates that the free-surfaces follows the velocity of the fluid, while 
Equation (10b) is a dynamic Dirichlet boundary condition for the Poisson prob- 
lem (10a) at each time t, whose initial condition is given by Equation (10c). 
A detailed explanation of the role of the term /j,(j) n in Equation (10b) is given 
in Appendix B. 

Introducing a non-physical reference domain (the Arbitrary Lagrangian Eu- 
lerian domain fi) and an arbitrary deformation p(x, t) such that p(Cl, t) = f2(t) 
as explained in Appendix A, Equations (10b) and (lOf) can be rewritten in 
terms of ALE derivatives 5/ St := d/dt + w ■ V as 

S ^ + (v-w)-V^=-gr ] + a s -x + ^\V(b\ 2 -^ n onT w (t) (11a) 

w ■ n = v ■ n on r(t), (lib) 

where w(p(x,t),t) = dp(x,t)/dt is the velocity of the domain deformation 
expressed in Eulerian coordinates (see Appendix A for the details). 

We remark that, whenever w — v, we recover a fully Lagrangian formulation: 
the ALE motion is, in this case, following the particles motion and 5 /St = D/Dt. 
Similarly, if the domain is fixed, and wc set w = 0, we recover the classical 
Eulerian formulation (10) on fixed domains, and 6/ St = d/dt. 

If the free-surface can be seen as the graph of a single valued function 
r](x,y,t) of the horizontal components x and y, then 

x = (x,y,r](x,y,t)) onr(t), (12) 

and, for a material particle on the free-surface, we have 

p(x,t)-e z =r](p(x,t),t), (13) 

where wc define rj(x,t) := r/(x ■ e Xl x ■ e y ,t). Taking the time derivative of 
Equation (13) we get 

^(x,t)-e z = ^(p(x,t),t) + ^(x,t)-V V (p(x,t),t) onf™, (14) 



that is, in Eulcrian form: 

Vz = ^L +v .Vr ] =^ onr(t), (15) 

where dr/fdz = 0. Proceeding in the same way for the ALE deformation and 
the ALE velocity of the domain, we get 

Wz = ^L +w .Vr ] = S £ onr(i). (16) 

Isolating drj/dt in Equation (15) and substituting in Equation (16), we get an 
alternative expression for Condition (lib) 1 , 

j- + (v - w) ■ Vr] = v ■ e z onTft), (17) 

which is valid for nonbreaking waves in which rj(x,y,t) is single- valued. 

Equation (17) is the kinematic boundary condition for the evolution of the 
unknown free-surface elevation rj{x, y, t) which is often found in the literature 
of semi-Lagrangian methods for potential free-surface flows [4, 5, 22]. 

Equation (17) is rather general, and is valid for arbitrary values of horizontal 
ALE velocities. Suitable values of V s and Vf can be plugged into the velocity 
field v = ( V s — Vf + V</>) to specify them for the desired reference frame. For 

instance, setting V s = Vf = 0, and w = ( 0, 0, S ) , one obtains 

^L = ^ + Vr? • (w - V6) := ^_^^Z_^^Z (is) 

St dz dz dx dx dy dy 

5 f t - -gV+l\V4>\ 2 + V<b-(w-V4) = -gr 1 +^- 1 -\Vcf>\ 2 , (19) 

which are the semi-Lagrangian equations written in an earth fixed reference 
frame, and null stream velocity used in [18]. 

In this work, we choose instead to solve the problem in a coordinate system 
attached to the ship hull. We move the reference frame according to the hori- 
zontal average velocity of the boat, that is we set Vf = (Vb • e x )e x =: — V^ = 
(— Voo,0, 0) and we assume V s = 0. With these values, Equations (17) and 
(11a) take the form 



Sij d<p 

Jt ~ dz 



Vrj-iw-V^-V^) (20) 



^ - -gv + llV^ + V^-iw-V^-V^-u^, (21) 



Conditions (lib) and (17) are equivalent in this framework. This comes from the obser- 
vation that, when T w (t) is the graph of the function r)(x,y,t), then the normal to the surface 
r™ (t) itself is given by a vector proportional to e z — V77. Substituting in Equation (lib) gives 
immediately Equation (17). 



which coincide with the ones proposed by Beck et al. [5], and in which the 
points on the free-surface move with an a priori arbitrary horizontal speed in 
the boat reference frame. 

2.2. Boundary integral formulation 

While Equation (10a) is time-dependent and defined in the entire domain 
il(t), we are really only interested in its solution on the boundary T(t). 

Knowledge of the solution on the free-surface part of the boundary is used 
to advance the shape of the domain in time, while Bernoulli's equation (3a) is 
used to recover the pressure distribution on the ship hull by postprocessing the 
solution on T h (t). 

In other words, at any given time instant t we want to compute <fi satisfying 



-A<f> = in n(t) (22a) 

4> = ~4> on T w (t) (22b) 

4> n =~<K on T h (t)yjT b (i)yjT f f(i) (22c) 

where 4> is the potential on the free-surface at the time t, and <f> n is equal to 
zero on T b (t) U T** (t) and to (V h (t) - V TO (t)) • n on T h (t). 

This is a purely spatial boundary value problem, in which time appears only 
through boundary conditions and through the shape of the time dependent 
domain. 

The solution of this boundary value problem allows for the computation of 
the full potential gradient on the boundary T(t), which is what is required in 
the dynamic and kinematic boundary conditions to advance the time evolution 
of both (f> and r\. 

Using the second Green identity 

f f Ou f f dv 

/ (-Au)vdx+ / — vds= / (-Av)udx+ / 7^-wds, (23) 

Jn J dn on J n J dn dn 

a solution to system (22) can be expressed in terms of a boundary integral 
representation only, via convolutions with fundamental solutions of the Laplace 
problem. 

We call G the fundamental solution, i.e., the function 

which is the distributional solution of 

- AG(x - x ) = S(x ) in R 3 



(25) 
where S(xq) is the Dirac delta distribution centered in a^- 



lim G(x — Xq) = 0, 

\x\— >oo 



If we select Xq to be inside il(i), use the defining property of the Dirac delta 
and the second Green identity, we obtain 

4>{x , t) = f [- (AG(x - a*)) <f>{x, t)] dQ = 

Jn(t) 

[ [(V<t>(x, t) ■ n)G{x - x ) - (VG(x - x ) ■ n)(j>{x, t)] AT. (26) 
Jv(t) 

In the limit for x touching the boundary T(t), the integral on the right 
hand side will have a singular argument, and should be evualated according to 
the Cauchy principal value. This process yields the so called Boundary Integral 
Equation (BIE) 



a<j> 

'r(t) 



on 



dr on T(t), (27) 



where a(x,t) is the fraction of solid angle Air with which the domain fi(£) is 
seen from x and the gradient of the fundamental solution is given by 

VG(r)-n = -—,. (28) 

4-7T \r\ 

The function a(x, t) can be computed by noting that the constant function 1 

is a solution to the Laplace equation with zero normal derivative, and therefore 

it must be 

f f)G 
a = - —dr on Fit), (29) 

Jr(t) dn 

in the Cauchy principal value sense. 

With Equation (27), the continuity equation has been reformulated as a 
boundary integral equation of mixed type defined on the moving boundary T(t), 
where the main ingredients arc the perturbation potential (j>(x, t) and its normal 
derivative (f> n (x,t). 

The domain deformation p{x,t) on the free-surface takes the form 

p(x,t) = (x,y,v(x,y,t)) onT w (t), (30) 

and one has to solve an additional boundary value problem to uniquely deter- 
mine the full ALE motion p(x,t). 

2.3. Arbitrary Lagrangian Eulerian motion 

When the arbitrary Lagrangian Eulerian formulation is used in the finite 
element framework (see, for example, [10]), the restriction to the boundary 
r(i) of the deformation p(x,t) is either known, or entirely determined by the 
equations of motion. In this case, an additional boundary value problem needs 
to be solved to recover the domain deformation in the interior of tt(t) starting 
from the Dirichlet values on the boundary. 
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Our situation is slightly different, since only the normal component of the 
motion is given on the boundary T(i), and we are not really interested in finding 
a domain motion in the interior of Q(t). 

In the dynamic and kinematic boundary conditions (20) and (21) we have 
the freedom to choose an ALE motion arbitrarily, as long as the shape of T(t) 
is preserved. In analogy to what is done in the finite element framework, we 
construct an additional boundary value problem to determine uniquely p(x,t). 

A typical choice in the finite element framework is based on linear elastic- 
ity theory, and requires the solution of an additional Laplace problem on the 
coordinates p(x,t), or, in some cases, a bi-Laplacian. This procedure can be 
generalized to surfaces embedded in three-dimensional space via the Laplace- 
Bcltrami operator. 

The Laplace-Beltrami operator can be constructed from the surface gradient 
'V s g(x,t), defined as 

V s a(x,t) :— Va — (Vo • n)n, Va s.t. a = a(x, t) on T(i), (31) 

where a is an arbitrary smooth extension oia(x, t) on a tubular neighborhood of 
T(i). Definition (31) is independent on the extension used (see, for example, [9]). 
Similarly, we indicate with V s the surface gradient computed in the reference 
domain T, with the same definition as in (31), but replacing x with x, and 
performing all differential operators in terms of the independent variable x 
instead of x. 

If we indicate with V s - the surface divergence (i.e., the trace of the surface 
gradient V s ), then the surface Laplacian A s and A s on T(t) and on T, are given 
by A s := V s • V s , and by A s := V s • V s . 

We use the shorthand notation ^""^{t) to indicate the intersection between 
T a (t) and T b (t), that is, 



7 «. 6 (t) = r«(*)nr»(*) a^b, (32) 

where a, b are either w, h, b or //. We indicate with 7(f) the union of all curves 
7 Q ' b (t). 

The curve j w,h (t) is usually referred to as the waterline on the hull of the 
ship. On j w - h (t), the domain velocity w has to satisfy the kinematic boundary 
condition for both the free-surface and the ship hull: 

w ■ n w = v ■ n w on 7 " , ' /l (t) 

h ( 33 ) 

wn h = on-r- h {t), 

where n w is the normal to the free-surface and n h is the normal to the hull 
surface. When both conditions are enforced, w is still allowed to be arbitrary 
along the direction tangent to the waterline. 

There are several options to select the tangent velocity w t defined as 

W t := w ■ (n h xn w ) = w t. (34) 
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A natural possibility is to choose zero tangential velocity. Other choices are 
certainly possible, and may be preferable, for example, if one would like to 
cluster computational nodes in regions where the curvature of the waterline is 
higher. In the experiments we present, the tangential velocity is always set to 
zero. 

Conditions (33) and zero tangential velocity, uniquely determine an evolution 
equation for the ALE deformation on j w - h . Here we summarize all boundary 
conditions for the evolution of p(x, t) on the entire 7: 



w ■ n h = 
w ■ (n h xn w ) = 



on j w - h (t) 

on -r- h (t) 

on -f w - h {t) 



(35a) 



w ■ n = v ■ n 

w ■ n^ = 
w ■ (n w x n ff ) = 



r Jf ( 

on j wJf (t) 



on 7""' J ■> (t) 



on 7 



(35b) 







on 7 



bjf 



(t). 



(35c) 



Expressing the boundary conditions (35) as a given velocity term w g , the evo- 
lution equation of j(t) become 



-^(x,t)=w g (p^(x,t)) on 7 

p 7 (x,0) =p (x) on 7. 



(36) 



A reconstruction of a reasonable ALE deformation on the entire T(t) is then 
possible by solving an additional elliptic boundary value problem, coupled with 
a projection on the surface of the ship hull and on the free-surface. Given a free- 
surface configuration r\ and the deformation p on the wireframe 7, in order to 
find a compatible deformation p on the entire T, we solve the additional problem 






-Ink 



onT 
on 7 



p = V h g 

p = 'Pr 1 g-=g + (v(g) - g ■ e z )e z 



onf 
onf u 



(37) 



p = g 



on f b U f // , 



where k{x) is the mean curvature of the domain L, i.e., the mean curvature of 
the hull on r and zero eveywhere else, while Vh is a projection operator on the 
hull surface. Similarly, V n is a (vertical) projection operator on the free-surface. 
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The auxiliary function g represents a surface that follows the waterline de- 
formation p . On the ship hull, it is a perturbation of the shape of the hull 
while everywhere else it is a minimal surface with boundary conditions imposed 
by p . On the free-surface, only its x and y components are used to determine 
p, while n (which satisfies the kinematic boundary conditions (17)) imposes the 
z component. 

2.4- Integro- differential formulation 

Putting everything together, the final intcgro-diffcrcntial system is given by 
the following problem: 

Given initial conditions 4>o & n d % on T u '(0), and p on 7, for each time 
t <G [0, T], find p, <j>, <j> n that satisfy 

r(t) o n Jv(t) on 

— = V<p((l),4> n ,n,w) 

Sii 

— = V v (<j>,(f> n ,r},w) 

<j>(x,Q) = <j>o(x) 
r)(x,0) = r] (x) 

<Pn = 4>n 

w ■ n = 

— A s g = —2nk 

dg 

-gj-(x,t) =w g (g(x,t)) 

p{g,o) =Po(9) 

p = Vg 



<P n G dr on T(t) 

(*) 


(38a) 


onTft) 


(38b) 


onTft) 


(38c) 


on T w (0) 


(38d) 


on r™(o) 


(38e) 


on T N {t) 


(38f) 


on r N (t) 


(38g) 


on r \ 7 


(38h) 


on 7 


(38i) 


on 7 


(38j) 


onf. 


(38k) 



13 



where we used the shorthand notations 



w : 


St 


(39a) 


v : 


= V 00 - Vf + V0 


(39b) 


r 4>(<P,<Pn,V, w ) ■ 


= (w - v) ■ V<f> - grj + aoo ■ x + -|V</>| 2 - n<f>„ 


(39c) 


v(<P,<Pn,V, w ) ■ 


= (w — v) ■ V77 + v ■ e z 


(39d) 


T N (t) : 


= r h (t)ur b (t)ur ff (t) 


(39e) 


4>n ■ 


_f(V b -Vf)-n on T h (t) 

\o onr b (t)urff(t), 

\Vh~9 onf' 1 


(39f) 


Vg: 


= \v il g:=g + {T 1 {g,t)-g-e z )e z onf 1 " 

\q onf fc uf^', 


(39g) 



and both the potential and the pressure in the entire domain can be obtained 
by postprocessing the solution to Problem (38) with the boundary integral rep- 
resentation (27) and with Bernoulli's Equation (3a). 

The full gradient of the perturbation potential on the surface Y{t) that ap- 
pears in Equations (39b), (39c) and (39d) is constructed from the surface gra- 
dient of 4> and from the normal gradient cj> n as 

V(p{x,t):=V s 4>{x,t)+<t> n {x,t)n. (40) 

A numerical discretization of the continuous Problem (38) is done on the 
fixed boundary Y of the reference domain f2, with independent variable x which 
will label node locations in a reference computational grid, and the motion 
p(x,t) will denote the trajectory of the computational nodes. 

3. Numerical discretization 

To approximate the continuous problem, we introduce a decomposition Yh 
of the domain boundary Y. Such partition is composed of three-dimensional 
quadrilateral cells (indicated with the index K), which satisfy the following 
regularity assumptions: 

1. Y = U{KeY h }; 

2. Any two cells K, K' only intersect in common faces, edges, or vertices; 

3. The decompositions Yh matches the decomposition Y = r w ur' l ur 6 ur^-' . 
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On the decomposition I 1 /,, we look for solutions (p h ,<j>h,(f>nh) in the finite di- 
mensional spaces Yh, Vh, and Qh defined as 

Y h := [u h G C°(f,0 3 | u fc |x G Py^) 3 , K G f fc } = span{«' fc }^ (41) 

^ := {^ G C°(f h ) | ^ G 7VW, K G f,,} = span{^}S (42) 

Qfc := {7ft G C°(f fc ) | 7h | X G P Q (if), tf G f ,,} = spanK}^. (43) 

Here, C°(r\) d is the space of continuous functions of d components over f\. 
Vy{K), Vv{K) and Vq(K) indicate the polynomial spaces of degree ry, ry 
and tq respectively, defined on the cells K. Finally, Ny, Ny and Nq denote 
the dimensions of each finite dimensional space. 

The most common approach for the solution of the boundary integral equa- 
tion (27) in the engineering community is the collocation boundary element 
method, in which the continuous functions <f> and Vc^-n are replaced by their dis- 
crete counterparts and the boundary integral equation is imposed at a suitable 
number of points on T(t). 

Once a geometric representation of the reference domain Th is available, we 
could in principle employ arbitrary and independent discretizations for each 
of the functional spaces Vh, Qh an d Y h . We choose instead to adopt an iso- 
parametric representation, in which the same parametrization is used to describe 
both the surface geometry and the physical variables. Thus, the deformed sur- 
face Th(t) (i.e., the map p(x,t)), the surface potential (f>h(x,t) and its normal 
derivative (f> n (x,t), are discretized on the fc-th panel as 

Qh = V h , Y h = V,f = span{<^e x , tp\ey, fUz}^., ( 44 ) 

where e x ,y,z are unit basis vectors identifying the directions of the x, y and 
z axis. We indicate with the notation {</>} and {4>n} the column vectors of 
time-dependent coefficients cf) 1 (t) and <p n J (t) such that 



N v 



h h {x,t) :=Y,<P i {t)'P i h{p- h 1 {x,t)) on T h (t) 



1 = 1 

N v 



(45) 



lh {x,t) :=Y,^{tW h {p- h \x,t)) on T h (t). 

3=1 



The map p h (x, t) is the inverse of the discretizaed ALE deformation, which 

reads 

N v 

p h (x,t):=J2x k (t)<p k h (x) onf, (46) 

fc=0 

where x k represents the time-dependent location of the vertices or control points 
that define the current configuration of T(t). 

To distinguish matrices from column vectors, we will indicate matrices with 
the bracket notation, e.g.., [M]. The ALE derivatives of the finite dimensional 
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<f>h(x,t) and p h (x,t) are time parametrized finite dimensional vectors in Vh and 
Yh, reading 

5 ^t)=j:%{tU{p- h \x,t)) 

l ;l , ( 47 ) 

w h (x,t) := S -^(x,t) = ]T ^(tM(pZ\x,t)). 

i—l 

Here, each function is identified by the coefficients {(/)}' and by the control points 
{a;}', where the ' denotes derivation in time. 

Finally, we can reconstruct the full discrete gradient V</>/j on Th(t) using the 
discrete version of the surface gradient V s and the normal gradient <p n h'- 

N v 

V<t>h(x, t) = V s 4>h{x, t) + (f> nh (x, t)n on T h (i). 

3.1. Collocation boundary element method 

The discrete version of the BIE, written for an arbitrary point y on the 
domain boundary, reads 



a(y,t)<j)(y,t) = 

m n v gG 

-EE^*) /. -Q-(y- xk ( u ^^))^l(u,v)J k (u,v,t)dudv 
fc=i i=i Jk n 

+ EE |W lG(y-x k (u,v,t)M(u,v)J k (u,v,t)dudv. (49) 
fe=i i=i ^ n ' J^ 

Here, we made use of the iso-parametric representation described in Appendix 
C, to decompose the integrals into the local contributions of the Nl basis func- 
tions in each of the M panels of the triangulation. 

The numerical evaluation of the panel integrals appearing in equation (49) 
needs some special treatment, due to the presence of the singular kernels G(y — 
x) and §^(y — x). Whenever y is not a node of the integration panel, the 
integral argument is not singular, and standard Gauss quadrature formulas are 
used. If y is a node of the integration panel, the integral kernel is singular and 
special quadrature rules are used, which remove the singularity by performing 
an additional change of variables (see, for example, [19]). In the framework 
of collocated BEM, an alternative possibility would have been represented by 
dcsingularizcd methods, in which the fundamental solutions are centered at 
points that are different from the evaluation points. Typically, this is obtained 
by centering the fundamental solutions at points that are slightly outside the 
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domain. Although these methods avoid dealing with singular integrals, they 
pose problems on establishing a general rule for suitable positioning of the fun- 
damental solutions centers. In the case at hand, the domain presents sharp 
edges and narrow corners (typically found at the bow or stern of a hull) which 
make the latter task nontrivial. 

Writing a boundary integral equation for each node x l , i — 1, . . . , Ny of the 
computational grid, we finally obtain the system 

[a]{<f>} + [N\{<t>} = [D]{<t>n}, (50) 

where 

• {4>} and {4> n } are the vectors containing the potential and its normal 
derivative node values, respectively; 

• [a] is a diagonal matrix composed by the a(xi(t)) coefficients; 

• [D] and [N] are the Dirichlet and Neumann matrices respectively whose 
elements are 

M . 

Dij = y2 / G(xi(t) -x k (u,v,t))ip J k (u,v)J k (u,v,t)dudv (51) 
k=i J& 

M f dG 

N i:j = Y^ L -Qn (x l (t)-x{u 1 v 1 t))ipl(u,v)J k (u,v,t)dudv. (52) 

(53) 

The evaluation of the nodal values for the solid angle fractions a, appearing 
in the BIE equation is obtained considering the solution to Laplace problem (22) 
when (j> = 1 in O(i). In such case, system (50) reads 

[a] {1} + [N] {1} = 0, (54) 

which implies 

N 

Oi = -J2 N iJ i = l,.-.,N. (55) 

i=i 

This technique is usually referred to as the Rigid Mode Technique, or Rigid 
Body Mode (RBM) (see, for example, [8]). It can be interpreted as an indirect 
regularization of the Neumann matrix [N] , which ensures that the matrix [a] + 
[N] possesses a zeroth eigenvalue corresponding to the rigid mode of the system. 

It is worth noting here, that there are arguments in favor of selecting discon- 
tinuous discrete representations for <p n = V(/> • n. Such considerations are based 
on the fact that surface normals are not necessarily continuous across neighbor- 
ing panels. In most cases though, the nodal {4> n } values obtained through the 
solution of system (50) represent a very reasonable approximation of the contin- 
uous normal potential gradient. But whenever the domain boundary presents 
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sharp features, such as corners or edges (see Fig. 2), a continuous approximation 
of {4> n } given by system (50) becomes extremely inaccurate, the exact normal 
potential gradient itself being not continuous. 




Figure 2: The two different normal unit vectors of a node placed on an edge 

To overcome such problem, in this work we employ a technique for the 
treatment of edges and corners, which was first developed by Grilli and Svendsen 
([14]). In this framework, the computational grid is first prepared so that on 
each edge separating two different boundary condition zones, the mesh nodes are 
duplicated. Hence, two distinct nodes are present in the mesh, each belonging 
to one of the two elements on the edge, and each having a different normal unit 
vector. In this way, the number of degrees of freedom of the system has been 
increased so as to allow discontinuous edge values for the normal derivative of 
the potential. 

3.2. SUPG stabilization 

For the time evolution of the dynamic and kinematic boundary conditions (38b) 
and (38c) , we need to evaluate the surface gradient of the veocity potential, as 
it appears in Equation (48). 

The gradient of the perturbation potential is not, in general, continuous 
across the edges of the panels that compose Th(t). As a consequence, the right 
hand side of equations (38b) is not single valued at the location of the computa- 
tional grid nodes. Thus, it is not possible to write a direct evolution equation for 
the nodes of the computational grid and for the corresponding potential values. 

A possible solution to this problem would be collocating the time evolution 
equations at marker points placed on the internal surface of each cell ([18]). 
Although this strategy would result in a single valued right hand side of equa- 
tions (38b), it would require that a new computational mesh is reconstructed 
from the updated markers position at each time step. A different approach 
consists in using using smooth finite dimensional spaces, as in [13]. This choice 
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would allow for the collocation of the evolution equations at the mesh nodes, 
but would in turn require the use of high order panels. 

In this work, wc choose instead to impose the evolutionary boundary con- 
ditions in weak form. Specifically, we employ an L 2 projection in the V/, space, 
to evaluate the right hand side of equations (38b) and (38c), namely 



St 



<p) =(Y*,<p) w v^eVk 



(|^) =0^)™ y ^ eVh > 



(56a) 
(56b) 



where 



(a,b) w = ab dr 

(a,b) = I ab dr. 
Jrit) 



(57) 




Figure 3: An example of the sawtooth instability developing on the stern of the hull without 
stabilization. 

This approach leads to extremely accurate estimations of the forcing terms in 
the free-surface evolution equations (39c) and (39d), with a very small computa- 
tional cost. Unfortunately, the equations right hand sides both contain transport 
terms (respectively V ??•(«? — Vcf) — V^ + Vf) and V <£•(«; — V0 — Voo + Vf)), 
that become dominant whenever (Vf — V^) is very different from V(/>. This 
causes a sawtooth numerical instability which in most cases develops in prox- 
imity of the hull stern, with consequent blow up of the simulations (an example 
of such instability is given in Figure 3). 

A natural, inexpensive and consistent stabilization algorithm which is able 
toreduce the observed instabilities is the Streamwise Upwind Petrov-Galerkin 
(SUPG) scheme (see, for example, [16, 25]). The SUPG stabilization consists in 
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replacing the plain L 2 projection in system (56) with the weighted projection 
'8<t> 



where 



<p + d-V,<p) =(V l p,<p + d-V s tp) w V(f£V h (58a) 

■£,ip + d-V,<p] =(V, l .r + d-V tlY ") ir V;tii, (581,) 



d:=r[^^). (59) 



r is a positive stabilization parameter which involves a measure of the local 
length scale (i.e. the "element length") and the local Reynolds and Courant 
numbers. Element lengths and stabilization parameters were proposed for the 
SUPG formulation of incompressible and compressible flows in [16], and an 
in depth study of the stabilization properties for free boundary problems was 
presented in [25]. However, to the best of the authors' knowledge, this is the 
first time that such stabilization technique is applied directly to a boundary 
value problem defined on a curved surface. 

Expressing system 58 in matrix form, we get 

[M][P r »M' = [PHTO (60a) 

[M][Pr™}{vY = [Pr-]{Vr,}, (60b) 

where the vectors and — sparse — matrices elements are computed by 

M ij :=(tp j ,ip i + d-V s tp i ) 

M r 

= y2 tfii{u,v)( y Lp i (u 7 v) + d-'V s tp t (u,v))j k (u,v,t)dudv (61a) 

k=i J * 

P-:=l S - **'(') <^(t) (61b) 

I otherwise 

Vi:={V^ + d-V s ^) 

M 

= y~] J V <t> (u,v)(ip i (u,v)+d-V s if i (u,v))J k (u,v,t)dudv (61c) 

k=i J * 

= y~] J V v (u,v)(ip i (u,v) + d-V s <p l (u,v))j k (u,v,t)dudv. (61d) 
k=i Jk 

3.3. Semi-discrete smoothing operator 

The semi-discrete version of the smoothing problem (37) can be obtained 
with a finite element implementation of the scalar Laplacc-Bcltrami operator, 
and its application to the different components of p. 
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The weak form of the Laplace-Beltrami operator on f of a scalar function u 
in V with Dirichlct boundary conditions u g on 7 and forcing term / reads 



V s u, V s f 



(/> <P)f 



on 7. 



(62) 



Here, Vb denotes the space of functions tp m V such that their trace on 7 is 
zero (see, for example, [6] and the references therein for some details on the 
numerical implementation of the Laplace-Beltrami operator). 
The semi-discrete form of Equation (62) can be written as 



[K]{u} = {F}, 



(63) 



where 



JT':= V S ^,V S <^ 



M . 

y~] / V s (p{(u,v) ■'V s Lp t k (u,v)J k {u,v)dudv 



**■=&&) 



(64) 



M 



Y] / /fc(u,w)<y2fc( u ^)^fc(M,w)dMdi;. 



fc=l 

The discrete Laplace-Beltrami operator is solved for the auxiliary vector 
variable g, whose finite dimensional representation is given by {g}. The forcing 
terms in the system are given by the mean curvature along the normal of L. In 
this case we write 

[K]{g} = {fc}, (65) 



where 



\K\ 



[K] ' 
[K] 
[K] 



n. 



(66) 



{k} := -2k 
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and the full matrix form of Problem (38) finally reads 
[a] {</>} + [N] {0} - [D] {<« 

[M][p TW }{<j>y-[p rw ]{Vt} 

[M}[p^]{ v y-[p^}{v,} 

[Pr-]{<K0)} - [Pr-KM 
[i¥-]{»7(0)} - [Pr»]W 



= 


(67a) 


= 


(67b) 


= 


(67c) 


= 


(67d) 


= 


(67e) 


= 


(67f) 


= 


(67g) 


= 


(67h) 


= 0, 


(67i) 



[I - Pr-KM - [I - Pr-]{M 

[p 7 ]{ 9 }' - [P 7 ]W 

[K]{g}-{k} 
[P]{g}-{x} 

where V is a numerical implementation of the projection operator (39g). On 
the hull surface, this operator is easily implemented analytically for simple hull 
shapes, such as that of the Wigley hull. In a more general case, it is desirable to 
have an implementation of the projection operator which can directly interrogate 
the CAD files describing the hull surface. In this work, we present results 
obtained with the former projection. A CAD based projection, which makes 
use of the OpenCASCADE library [1], is currently being implemented. Some 
work is still required though, to render our full discretization robust with respect 
to arbitrary hull geometries. 

3.4- Time discretization 

System (67) can be recast in the following form 

F(t,y,y') = 0, (68) 

where we grouped the variables of the system in the vector y: 



y={ W )■ (69) 




Equation (68) represents a system of nonlinear differential algebraic equa- 
tions (DAE), which we solve using the IDA package of the SUNDIALS Open- 
Source library [15]. As stated in the package documentation (see p. 374 and 375 

in [15]): 2 

The integration method in IDA is variable-order, variable-coefficient 
BDF [backward difference formula], in fixed-lcading-cocfficient form. 



2 We quoted directly from the SUNDIALS documentation. However, we adjusted the no- 
tation so as to be consistent with ours and we numbered equations according to their order 
in this paper. 
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The method order ranges from 1 to 5, with the BDF of order q given 
by the multistep formula 

<? 
^a^iVn-i = h n y n , (70) 

where y n and y n are the computed approximations to y(t n ) and 
y'{t n ), respectively, and the step size is h n — t n — t n -\- The coeffi- 
cients a n .i are uniquely determined by the order q, and the history 
of the step sizes. The application of the BDF [in Eq. (70)] to the 
DAE system [in Eq. (68)] results in a nonlinear algebraic system to 
be solved at each step: 

R{Vn) = ■H*n)3/T»ftn i y^ a n,t3/n-i) =0. (71) 

Regardless of the method options, the solution of the nonlinear sys- 
tem [in Eq. (71)] is accomplished with some form of Newton itera- 
tion. This leads to a linear system for each Newton correction, of 
the form 

J[y n , m +i - y n , m ] = -R(y n ,m), (72) 

where y n . m is the mth approximation to y m . Here J is some approx- 
imation to the system Jacobian 

j = dR = OF + a dF 
dy dy dy' ' 

where a = a n fi/h n . The scalar a changes whenever the step size or 
method order changes. 

In our implementation, we assemble the residual R(y n ,m) at each Newton correc- 
tion, and let SUNDIALS compute an approximation of the Jacobian in Eq. (73). 
The final system is solved using a preconditioned GMRES iterative method (see, 
e.g., [12]). 

Despite the increase in computational cost due the implicit solution scheme, 
the DAE system approach presents several advantages with respect to explicit 
resolution techniques. First, it is worth pointing out that among the unknowns 
in equation (69), the coordinates of all the grid nodes (except for the vertical 
coordinates of free-surface nodes) appear in the DAE system as algebraic com- 
ponents. Their evolution is in fact not described by a differential equation, but 
computed through the smoothing operator. Yet, the time derivative of such co- 
ordinates, i.e.: the ALE velocity w is readily available through the evaluation 
of the BDF (70). Such velocites are used in the differential Equations (67b) and 
(67c) which appear in the DAE system. In particular, at each time step, the 
convergence of Newton corrections (72) ensures that the vertical velocity of the 
nodes is the one corresponding to the w velocity originated by the horizontal 
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nodes displacement computed by the smoothing operator. In a similar fashion, 
the DAE solution algorithm computes the ALE time derivative of the velocity 
potential at each Newton correction. Such derivative is plugged into Bernoulli's 
equation (3a) to evaluate the pressure on the whole domain boundary F(£), 
without requiring the solution of additional boundary value problems for -^. 
Finally, the resulting pressure field is integrated on the ship wet surface T h (t) 
to obtain the pressure force acting on the ship. Since this operation is done 
at the level of each Newton correction, possible rigid motions of a hull along 
its six degrees of freedom can be accounted for in a very natural way in the 
DAE framework, by adding the six differential equations of motion governing 
the unknown rigid displacements to the DAE system. The latter — strongly 
coupled — fluid structure interaction model is currently under development and 
results will be presented in future contributions. 

3.5. Adaptive mesh refinement 

There are two main advantages of using a Galerkin formulation for the evo- 
lution equation of the free-surface, as well as for the computation of the full 
ALE deformation on the surface T(t). On one hand, it allows the use of fully 
unstructured meshes, with an immediate simplification in the mesh generation. 
On the other hand, it enables the use of a wide set of local refinement strate- 
gies based on a posteriori error estimators, which are rather popular in the 
finite element community. As a result, the mesh generation and adaptation are 
fully automated, and the resulting computational grids are shaped based on the 
characteristic of the problem solution, rather than on a-priori heuristic choices. 

In this work, we use a modification of the gradient recovery error estimator 
by Kelly, Gago, Zienkiewicz and Babuska [17, 11], a choice mostly motivated by 
its simplicity (see [2] for more details on this and other error estimators). 




Figure 4: A mesh refinement step. 

At fixed intervals in time, the surface gradient of the finite element approxi- 
mation of cj> is post-processed. This provides a quantitative estimate of the cells 
in which the approximation error may be higher. In particular, for each cell K 
of our triangulation we compute the quantity 



'K 



h 
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<)K 



nax] dj, 



(74) 



where [V s </> • u^k] denotes the jump of the surface gradient of <fi across the edges 
of the triangulation element K. The vector tiqk is perpendicular to both the 
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cell normal n, and to the boundary of the element K, and h is the diamater of 
the cell itself. 

Roughly speaking, tk gives an estimate of how well the trial space is ap- 
proximating the surface gradient of (p. The higher these values are, the smaller 
the cells should become. The estimated error per cell tk is ordered, and a fixed 
fraction of the cells with the highest and lowest error tk are flagged for refine- 
ment and coarsening. The computational grid is then refined, ensuring that any 
two neighboring cells differ for at most one refinement level. 

Standard interpolation is used to transfer all finite dimensional solutions 
from one grid to another, and a geometrically consistent mesh modification al- 
gorithm is used to collocate the new nodes coordinates as smoothly as possible 
(see [6] for a detailed explanation of this algorithmic strategy). The result- 
ing computational grid is non conformal. At each hanging node, all the finite 
dimensional fields are constrained to be continuous. This results in a set of 
algebraic constraints to be applied to all the degrees of freedom associated with 
the hanging nodes, which are eliminated from the final system of equations via 
a matrix condensation technique. 

Most of these algorithmic strategies are based on the ones which were al- 
ready implemented in the deal. II finite element library [3] for trial spaces of 
finite elements defined in two and three dimensions, and were modified to allow 
their use on arbitrary surfaces embedded in three dimensions. An example of a 
refinement step is presented in Figure 4. 

After each coarsening and refinement step, the system of differential alge- 
braic equations is restarted with the newly interpolated solution as initial con- 
dition. A state diagram for the entire solution process is sketched in Figure 5. 

4. Numerical simulations and results 

The test case presented in this work is the problem of a Wigley hull ad- 
vancing in calm water at speed V^ parallel to the longitudinal axis, with fixed 
sinkage and trim. In naval engineering, the Wigley hull is commonly used as 
a benchmark for the validation of free-surface flow simulations. This is mainly 
due to its simple shape, defined by the equation 



y(x,z) = — 



ZN-' 

f 



(75) 



In our simulations the boat length, beam and draft values used are respectively 
L = 2.5 m, B = 0.25 m, and T = 0.15625 m. A sketch of the resulting hull 
shape is presented in Fig. 6, which represents a set of vertical sections of the 
hull. 

In the numerical simulation setup, the boat is started at rest, and the its 
velocity is increased up to the desired speed Voo with a linear ramp. The 
simulation is then carried on until the flow approaches a steady state solution. 
The presence of the ramp is not needed for the stability and convergence of the 
solution, which are also obtained imposing an impulsive start of the water past 
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Figure 5: State Diagram 

the hull. Still, inducing slower dynamics in the first seconds of the simulations, 
the linear ramp allows for higher time steps and faster convergence. 

To compare the non linear free-surface BEM solutions with the experimental 
results presented in [20], we considered six different surge velocities Vco, corre- 
sponding to the Froude numbers reported in Tabic 1. For each of these Froude 
numbers, numerical solutions were obtained using a refined and a coarse mesh, 
in which the adaptive mesh refinement algorithm was tuned in order to main- 
tain the cell dimensions over a given minimum value, and in order to limit the 
number of degrees of freedom under a maximum value. 

Despite the fact that the BEM solver developed allows for the use of panels 
with arbitrary order, the results we present only refer to iso-parametric bilinear 
elements. The proposed stabilization mechanism seems insufficient for higher 
order quadrilateral elements, especially when high Froude numbers are consid- 
ered. 

We are currently investigating alternative stabilization procedures as well as 
finer tuning of the SUPG parameters to increase the robustness of the algorithm 
when high order panels are used. 

A contour of the wave elevation field for the regime solution obtained at the 
various Froude numbers is presented, along with the final mesh, in Figures 7, 
8 and 9. The pictures show how the adaptive mesh refinement leads to an 
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Figure 6: Vertical sections of the Wigley hull used for the simulations, generated by planes 
normal to the longitudinal axis of the hull. 



Voo 


1.2381 ^ 

S 


1.3223 *2 
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1.4312 2i 

S 


1.5649 ^ 

S 


1.7531 ^ 

S 


2.0205 ^ 

S 


Fr 


0.250 


0.267 


0.289 


0.316 


0.354 


0.408 



Table 1: Wigley hull surge velocities imposed in each numerical simulation, and the corre- 
sponding Froude numbers Fr = SSL. . 



automatic clustering of mesh cells in regions with highest solution gradients. In 
this way the numerical solutions capture in a very accurate manner the physical 
characteristics of the wave patterns, requiring a very limited number of degrees 
of freedom. The biggest final mesh is in fact only composed by roughly 6000 
nodes, but it allows for a very good reconstruction of the Kelvin wake, extending 
for several wavelengths past the surging hull. 

The simulations required 12 hrs (for the coarse meshes) to 48 hrs (for the re- 
fined meshes) to reach the steady state solution, on single SMP nodes of the Arc- 
tur cluster of the Italian/Slovenian interstate cooperation Exact-Lab/Arctur. 

The wave profiles on the surface of the Wigley hull obtained with the present 
method, are compared with the corresponding experimental results in Fig. 10. 
In each plot, the abscissae represent the dimensionless coordinate x/L along the 
boat, while the ordinates are the dimensionless wave elevations rj' = t/%. For 
all the Froude numbers considered, the method presented seems able to predict 
qualitatively correct wave profiles. Moreover, the wave elevation in proximity of 
the bow of the boat is reproduced with very good accuracy in all the test cases 
considered. On the other hand, in all the numerical curves, we observe a small 
spatial oscillation superimposed to the main wave profile. The wave length 
of such oscillation seems proportional to the local mesh cells size, while the 
amplitude is slightly higher for finer meshes. This suggests that better tuning 
of the SUPG stabilization parameter may be needed for this kind of boundary 
value problems. 
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(b) Fr = 0.267 
Figure 7: Mesh refinements and contours (I). 



5. Conclusions 



An accurate and efficient boundary element method for the simulation of 
unsteady and fully nonlinear potential waves past surging ships was developed, 
implemented and tested. Compared to existing algorithms, the method presents 
several innovative features which try to address some of the most CPU intensive 
aspects of this kind of computations. 

The most innovative idea behind the proposed method is the fact that the 
equations are studied on a fixed reference domain, which is deformed through 
an arbitrary Lagrangian Eulerian map that keeps track of the physical shape of 
the water domain around the ship. Some aspects of this approach resemble the 
semi-Lagrangian formulation of the potential wave equations, but here they are 
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(a) Fr = 0.289 
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(b) Fr = 0.316 
Figure 8: Mesh refinements and contours (II). 

taklcd using powerful differential geometry tools, combined with finite clement 
techniques for arbitrary surfaces. 

This reformulation in terms of a fixed reference domain presents severe sta- 
bility issues in presence of a forward ship motion, or in presence of an incident 
stream velocity. Stabilization is achieved via a weighted SUPG projection, which 
allows the use of fully unstructured meshes, and guarantees an accurate recon- 
struction of the velocity fields on the mesh nodes, also when low order finite 
dimensional spaces are used for the numerical discretization. 

To the best of the authors' knowledge, such formulation has never been suc- 
cessfully used in ship hydrodynamic problems in presence of a non zero stream 
velocity, due to the free-surface instabilities. 

With respect to existing methods, the combination of the semi-Lagrangian 
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(a) Fr = 0.354 




(b) Fr = 0.408 
Figure 9: Mesh refinements and contours (III). 

approach with the SUPG stabilization eliminates the need for periodic remesh- 
ing of the computational domain, and opens up the possibility to exploit local 
adaptivity tools, typical of finite element discretizations. 

We exploit these ideas by employing simple a-posteriori error estimates to 
refine adaptively the computational mesh. Accurate results are obtained even 
when using a very limited number of degrees of freedom. 

Implicit BDF methods with variable order and variable step size are also 
employed, which render the final computational tool very attractive in terms of 
robustness. 

A direct interface with standard CAD file formats is currently under devel- 
opment, and our preliminary results indicate that the final tool could be used 
to study the unsteady interaction between arbitrary hull shapes and nonlinear 
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Figure 10: Comparison of predicted water profiles with the University of Tokyo experimental 
results (- -*- -). Both coarse mesh (— -) and refined mesh results ( — ) are shown in the plots. 



water waves in a robust and automated way. 
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Appendix A. Arbitrary Lagrangian Eulerian Formulation 



x 




Figure A. 11: Schematic representation of a Lagrangian motion and of an arbitrary Lagrangian 
Eulerian motion. 



A motion is a time parametrized family of invertible maps which associates 
to each point A in a reference domain £1 its position in space at time t: 



p:!lxl4 



(x,t) 1-4 x = p(x,t). 



(A.l) 



The domain tt(t) at the current time can be seen as the image under the 
motion p of a reference domain Cl, i.e., p(Ct, t) = il(t). We will indicate with the 
: symbol a material motion, or a motion in which the points x label material 
particles. 

If one does not want to follow material particles with the domain fl(t), it 
is possible to introduce another intermediate motion, called the ALE motion, 
with which we represent deformations of the domain Q(t): 



p:!lxl4 



(x,t) 1-4 x = p(x,t). 



(A.2) 



These motions can be rather arbitrary, as long as the shape of the domain 
Ci(t) is preserved by the motion itself, i.e., p(fl, t) = O(t). The points labelled 
with x do not, in general, represent material particles. See Figure A. 11 for a 
schematic representation of a Lagrangian motion and of an ALE motion. 
Given a Lagrangian held q : O x K i— >• R, its Eulerian representation is 



q(x,t):=q(p \x,t),t), Vx€fl(t), 



(A.3) 



:V2 



while, given an Eulerian field q : Cl(t) xl^I, its Lagrangian representation 
would be 

q(x, t) := q(p(x, t), t), \fx g 0. (A.4) 

A similar structure is valid for ALE fields: 

q(x,t):=q(p- 1 {x,t),t) 7 Vx g fi(t), (A.5) 

q(x, t) := q(p(x, t),t), Vx g f2. (A.6) 

The fluid particle velocity v which appears in Problem (1) is the Eulerian 
representation of the particles velocity 

v(p(x,t),t) = v(x,t):=^^-. (A.7) 

In a similar way, we define the Eulerian representation of the domain velocity, 
or ALE velocity the held w such that 

w(p(x,t),t)=w(x,t):= dp ^ t \ (A.8) 

Time variations of physical quantities associated with material particles are 
computed at hxed x, generating the usual material derivative 



Dq(x,t) dq(p(x,t),t) 



Dt dt 



--P-Hv.t) dt 



dq{x ' t] \-v-Vq(x,t). (A.9) 



In a similar fashion, if we compute quantities at fixed ALE point x, we 
obtain the ALE time derivative, wich we will denote as 



Sq(x,t) _ dq(p(x,t),t) 
St "~ dt 



= d ^^+wVq(x,t). (A.IO) 

x=p- 1 (x,t) 0t 



The ALE deformation allows one to define the equations of motions in Prob- 
lem (I) in terms of the fixed ALE reference domain 17, while still solving the 
same physical problem. On the free-surface part of the boundary, the ALE 
motion needs to follow the evolution of the fluid particles in order to maintain 
the correct shape of the domain Q(t), in particular the minimum requirement 
for the ALE motion on the free-surface is given by the free-surface kinematic 
boundary condition 

wn = vn onTft), (A.ll) 

which complements the dynamic boundary condition (lc), and provides an 
evolution equation for the normal part of the ALE motion on the boundary 
T w (t) . In terms of the ALE motion, the above condition reads 

-£(x,t)-n = v(p{x,t),t)-n onT w (t). (A.12) 
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Equations (lc) and (A. 11) are usually referred to as the free-surface dynamic 
and kinematic boundary conditions, since they represent the physical condition 
applied to the free-surface (equilibrium of the pressure on the water surface) and 
its evolution equation (the shape of the free-surface follows the velocity field of 
the flow). 

Appendix B. Numerical beach 

A draw back of using an homogeneous Neumann boundary conditions for 
the vertical far field boundary condition is that it reflects energy back in the 
computational domain. 

We use an absorbing beach technique (see, for example, [7]), in which we 
add an artificial damping region away from the boat, used to absorb the wave 
energy. A damping term can be seen as an additional pressure P acting on the 
free surface. In such case, Bernoulli equation on the free-surface becomes 

^ +gz -a s -x+ l -\V<S>\ 2 + -=0 on T w (t), (B.l) 

and one can show that the resulting rate of energy absorption is 

1r = / r . p *" dr - < R2 » 

A natural way to construct the damping pressure P is then to use a term which 
is proportional to the potential normal derivative </>„, which grants a positive 
energy absorption at all times. 

The dynamical free-surface boundary condition modified to account for the 
damping term reads 

^+gz-a s -x+^\V^>\ 2 -^ n = on T w (t), (B.3) 

where 

^ /ma^p)\ 2 _ (B4) 

and Xd is the x coordinate value in which the artificial damping starts to act, 
while Ld is the length of the absorbing beach, as in Figure 1. 



Appendix C. Isoparametric discretization 

We approximate the geometry of the domain boundary by means of arbitrary 
order quadrilateral panels. The edges of each panel are defined by polynomial 
functions, while their internal surface is described by polynomial tensor prod- 
ucts. 
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Figure C.12: A linear panel and the reference element. The red dots represent the additional 
degrees of freedom of a quadratic panel. 



In particular we employ the Lagrangian shape functions Ni(u,v) 1 = 1,..., N^ 
on the reference panel (Fig. C.12, on the left). Thus, the local parametrization 
of the k-th panel reads 



N, 



x k {u,v):^^x kl Ni{u,v) u,ve[0,l} 2 

x k {u,v,t) :^^2x kl {t)Ni{u,v) u,v€ [0,1] 2 . 



(C.l) 



1=1 



The parametrization weights in Eq. (C.l) are the positions of the nodes in the 
reference domain T^, or in the current domain Th(t). ki denotes the local to 
global numbering index which identifies the Nl basis functions ip kl which are 
different from zero on the A:-th panel. 

The global basis functions (p l (x) can be identified and evaluated on each 
panel K via their local parametrization as 



N L 



: ~Pk 



l k (u,v) :=<p l (x k (u,v)) = ^25 lkl Ni(u,v), 



Sn = 



1=1 



1 if i = j 

otherwise. 



(C.2) 



In this framework, the local representation of 4>{x k {u, v, t),t) and of its nor- 
mal derivative on the fc-th panel are 



N L 



N 



cf> k (u,v,t) = J2(f> k '(t)N l (u,v) (f>„ k (u,v,t) = J2<t> n kl (t)N l {u,v), (C.3) 



1=1 



1=1 
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where <p l ,<j) n l I = 1, . . . ,Nl are the nodal values of the potential and of its 
normal derivative in panel k. On each point of the panel, it is possible to 
compute two vectors tangential to T^(t) as 

%{u,v,t)=^2x k '{t)-^-(u,v) t k v (u,v,t)^J2 xkl ^^(^ v ^ ( C - 4 ) 
1=1 u 1=1 

from which the external normal unit vector n is obtained as 
n k vect (u,v,t)=t k u (u,v,t)xt k (u,v,t), n k (u,v,t)= ^^% - (C.5) 

The same can be done for vectors tangential and normal to T, by simply 
replacing x kl (t) with x ' in the definitions above. We will denote those vectors 

with i u (u,v), i v (u,v), n vect (u,v) and h (u,v), respectively. 

Integrals on a panel K (or K), can be computed in the reference domain 
[0, l] 2 , by observing that dr = J k (u, v, t) dudv, where J k (u,v,t) := \n k ect (u,v,t)\ 
(or dT = J k (u, v)dudv, where J k (u,v) :— \h vect (u,v)\). 

We finally introduce the following differential operators 

D k (u,v,t) := V uv x k {u,v,t) gR 3x2 

G k (u, v, t) := D k (u, v, t) T D k (u, v, t) E K 2x2 

D k (u, v) := V uv x k (u, v) e R 3x2 

G k (u, v) := D k (u, v) T D k (u, v) € M 2x2 , 



(C.6) 



where G k and G k are the first fundamental forms in the element K and K. 
Making use of such operators, we can express the local surface gradients of the 
global basis functions as 

V s tp\p~ l (x,t)) | x=Bfc ( u ,„, t ) = DkiGk^V uv Lp\(u,v,t) 

=:V a ip%{u,v,t) 

V s ip l {x) \x=X k (u,v) = I>fe(Gfe) _1 V ttt ,(y9 4 fc (u,w) 

=: V s ipl(u,v), 

which we will indicate with the same symbol as the spatial surface gradients. 
The surface gradient of a finite dimensional vector can then be computed by 
Equation (48). 

Appendix D. Nomenclature 
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